Generalized Function
- Function Purpose
- Simulation Workflow
  - Decision Options:
Estimates
Update based on field observations
Working with data from Bavi
- Extract relevant information from the data
- Updating with parts_used data
Bayesian Updating
- Check for species consistency
  - Species not included in first rounds
  - Species not included in later rounds
Generalize the update to all species in the collection (from all trips)

This project is concerned with the interrelated issues of conservation and nutrition. Specifically we are applying Monte Carlo simulations based on expert knowledge to assess the dietary and conservation implications of allowing forest access to communities around a conservation area in Vietnam.

We use well-established and flexible method for decision modeling with expert-derived probabilities and variances to estimate outcomes. With this method we can capture the variability and uncertainty in model inputs. We can apply this to guide decision-making by identifying likely outcomes and their variability.

We go through several iterative processes to simulate the effects of different park management strategies (e.g., restricting access, providing seedlings) on species conservation, dietary diversity, financial cost to the park service and cultural heritage, i.e. retaining or losing traditional knowledge.

Should Bavi Park Authority and Local Gov’t authorities be willing to allow access to the buffer zone or maintain current mostly exclusionary practice - Implications for sustainable food environments of local people

Generalized Function

The function simulates outcomes for decision options. The simulate_outcomes_for_decisions function is a Monte Carlo simulation framework designed to evaluate multiple decision options against predefined outcomes, incorporating stochastic variability (variability_list) to reflect real-world uncertainties. This function is intended for use in scenarios where decisions have probabilistic impacts on outcomes, such as ecological management, conservation planning, or policy evaluation. Source the simulate_outcomes_for_decisions function with source.

source("functions/simulate_outcomes_for_decisions.R")

Function Purpose

The function calculates the expected outcomes for a series of decision options under uncertainty. Each decision option is evaluated against a set of outcomes using a weighted linear model, with added random variability_list to account for variability. The outputs are stored in a structured format for further analysis.

Function Inputs include:

decision_options, a list where each element is a named vector representing a specific decision option. Each vector contains binary values (1 or 0), indicating whether a particular decision variable is active or inactive. For example: control_access = c(1, 0, 0, 0, 0) indicates that control_access is the first option in the list of decision options.

outcome_effects, a matrix where rows represent outcomes (e.g., conservation success, financial cost). The columns correspond to decision variables. The elements specify the weight or influence of a decision variable on a given outcome. i.e. matrix(c(0.8, 0.6, 0.5, 0.6, 0.7),nrow = 1) for eh species Conservation outcome_effects corresponding to the columns of decision_options.

The variability_list is a named list of functions, with each function generating random noise specific to a specific outcome. This allows for custom noise distributions to reflect expert knowledge or empirical data on variability. i.e. list(SpeciesConservation = function(n) rnorm(n, mean = 0, sd = 0.05)) is a function assigning the noise to the interaction between species conservation and the decision options.

n_simulations is the number of Monte Carlo iterations to perform for each combination of decision and outcome.

The function returns results as a nested list with: Decision options (e.g., control_access). Outcome names (e.g., SpeciesConservation). Values: A vector of simulated outcome values for each combination of decision and outcome.

Simulation Workflow

Decision Options:

The function iterates over each decision option in decision_options. For each decision option, the function calculates all outcomes defined in the outcome_effects matrix. Outcomes are calculated as a weighted sum of decision variables (defined by the outcome_effects matrix) plus a random variability_list term. variability_list is generated, allowing for outcome-specific variability. Results are stored in a nested list, enabling easy access for further analysis.

For a conservation project, this function will be used to simulate the effects of different park management strategies (e.g., restricting access, providing seedlings) on:

Species conservation: Probability of preserving species.
Dietary diversity: Impact on local food systems.
Expenditure: Financial cost to the park service.
Cultural heritage: Loss of traditional knowledge.

By comparing the simulated outcomes for each strategy, decision-makers can identify trade-offs and prioritize actions. The function can handle multiple decision options and outcomes. It incorporates expert-defined variability_list distributions for outcome variability. It also explicitly separates decision logic, outcome effects, and variability_list, ensuring clarity and reproducibility. For each decision option and outcome, the function calculates the Outcome Value using the following formula:

\[ \text{Outcome Value} = \sum_{i=1}^{n} \text{Decision}_i \times \text{Weight}_i + \text{variability_list} \]

Where

$\text{Decision}_i$: A binary value ($1$ or $0$) indicating if the $i$-th decision is active.
$\text{Weight}_i$: The influence of the $i$-th decision variable on the outcome.
$\text{variability_list}$: A random term representing variability or uncertainty, specific to the outcome and decision, as defined by the variability_list.

# Load necessary libraries
library(DiagrammeR)
# library(DiagrammeRsvg)
library(rsvg)

# Create the updated graph
forest_access_graph <- grViz("
  digraph forest_access_model {
    graph [layout = dot, rankdir = LR, splines = line]

    // Define node styles
    node [shape = box, fontname = Helvetica]

    // Define nodes
    ControlAccess [label = 'Control Access']
    ProvideAccessControlHarvest [label = 'Provide Access + Control Harvest']
    AllowAccessForestParts [label = 'Allow Access to Forest Parts']
    AllowAccessLimitedTime [label = 'Allow Access for Limited Time']
    ProvideSeedlings [label = 'Provide Seedlings']

    SpeciesConservation [label = 'Species Conservation']
    DietaryDiversity [label = 'Dietary Diversity']
    SustainableFoodEnv [label = 'Sustainable Food Environment']
    ParkExpenditure [label = 'Park Expenditure']
    TraditionalKnowledge [label = 'Traditional Knowledge']

    // Define relationships between nodes
    ControlAccess -> {SpeciesConservation DietaryDiversity SustainableFoodEnv ParkExpenditure TraditionalKnowledge}
    ProvideAccessControlHarvest -> {SpeciesConservation DietaryDiversity SustainableFoodEnv ParkExpenditure TraditionalKnowledge}
    AllowAccessForestParts -> {SpeciesConservation DietaryDiversity SustainableFoodEnv ParkExpenditure TraditionalKnowledge}
    AllowAccessLimitedTime -> {SpeciesConservation DietaryDiversity SustainableFoodEnv ParkExpenditure TraditionalKnowledge}
    ProvideSeedlings -> {SpeciesConservation DietaryDiversity SustainableFoodEnv ParkExpenditure TraditionalKnowledge}
  }
")

forest_access_graph

In this model ParkExpenditure refers to the additional costs incurred (or saved) by implementing food access strategies, more is worse. All other values are positive, more is better.

Estimates

Add species with simplified initial input

Use the add_species_data.R function with other functions from the tidyverse

source("functions/add_species_data.R")
library(tidyverse)

The Effect represents the magnitude and direction of influence that a specific decision option has on an outcome for a given species. Positive values indicate a beneficial effect and negative values visa versa. For instance, a value of 0.8 for control_access on Species Conservation implies a strong positive effect, whereas -0.2 for provide_seedlings on Dietary Diversity indicates a slight reduction in dietary diversity. A value of 0 means no effect. The Effect value represents the proportional change in the outcome relative to its current state for a given species. The measure allows for a direct comparison of how different interventions influence species-level outcomes.

The SD (also Variability_SD) quantifies the uncertainty or variability in the effect values. A smaller SD (e.g., 0.05) indicates high confidence in Effect, suggesting that the impact is consistent across scenarios. Larger SD (e.g., 0.2) reflects greater uncertainty, for example, an SD of 0.05 for control_access on Species Conservation indicates low variability and high confidence, whereas an SD of 0.15 for provide_seedlings on Dietary Diversity suggests moderate uncertainty.

Decision options for sustainable forest food

Define decision_options for managing forest access (limited time, area or volume of harvest) or species access (seedlings provided)

# Define decision options and outcomes
decision_options <- c("control_access", "provide_access_control_harvest", 
                      "allow_access_forest_parts", "allow_access_limited_time", "provide_seedlings")

# Define readable labels
decision_labels <- c(
  "control_access" = "Control Access",
  "provide_access_control_harvest" = "Regulated Harvesting",
  "allow_access_forest_parts" = "Forest Access",
  "allow_access_limited_time" = "Limited Access",
  "provide_seedlings" = "Seedling Provision"
)

Outcomes of sustainable forest food decisions

Define decision outcomes for diets and conservation

outcomes <- c("SpeciesConservation", "DietaryDiversity", "SustainableFoodEnv", 
              "ParkExpenditure", "TraditionalKnowledge")

outcome_labels <- c(
  "SpeciesConservation" = "Species Conservation",
  "DietaryDiversity" = "Dietary Diversity",
  "SustainableFoodEnv" = "Sustainable Food Environment",
  "ParkExpenditure" = "Park Expenditure",
  "TraditionalKnowledge" = "Traditional Knowledge"
)

Species Classification

Below is the list of species mentioned and observed as useful in the field. The list is organized by functional role of species in conservation, food security. Although we are mainly concerned about forest conservation in protected forests we also include information about the role of the species in agroforestry systems since many of the species are kept in local homegardens by those community members who have access and land.

Mushrooms:
- Auricularia heimuer (Jelly ear fungus)
- Volvariella volvacea (Paddy straw mushroom)
Grasses & Bamboo:
- Ampelocalamus patellaris (Bamboo)
- Arundinaria baviensis (Bamboo)
- Bambusa blumeana (Thorny bamboo)
- Coix lacryma-jobi (Job’s tears, grain)
- Dendrocalamus velutinus (Bamboo)
Ferns:
- Diplazium esculentum (Edible fern)
- Nephrolepis cordifolia (Sword fern)
Shrubs & Small Trees:
- Bischofia javanica (Bishop wood, tree)
- Canarium bengalense (Canarium tree)
- Canarium nigrum (Black canarium tree)
- Canarium tramvangum (Canarium tree)
- Euodia bodinieri (Rutaceae tree/shrub)
- Maesa perlarius (Shrub)
- Schefflera heptaphylla (Umbrella tree)
Climbers & Vines:
- Cissus hastata (Vine)
- Erythropalum scandens (Climber)
- Hodgsonia macrocarpa (Climber, liana)
- Paederia lanuginosa (Climber)
- Passiflora foetida (Passionflower vine)
Herbs & Vegetable Crops:
- Alpinia xanthioides (Ginger family)
- Begonia aptera (Begonia, herb)
- Begonia tonkinensis (Begonia, herb)
- Bidens pilosa (Blackjack, medicinal herb)
- Centella asiatica (Gotu kola, medicinal herb)
- Costus speciosus (Spiral ginger, herb)
- Crassocephalum crepidioides (Edible herb)
- Curculigo orchioides (Orchioides, herb)
- Dioscorea hamiltonii (Yam, tuber crop)
- Gynura procumbens (Longevity spinach, herb)
- Gynura sarmentosa (Herb)
- Melientha suavis (Pak wan, vegetable crop)
- Momordica cochinchinensis (Gac fruit, climbing vegetable)
- Peperomia pellucida (Shiny bush, herb)
- Physalis angulata (Ground cherry, herbaceous plant)
- Piper lolot (Wild betel leaf, herbaceous vine)
- Plantago major (Plantain herb)
Trees:
- Baccaurea ramiflora (Burmese grape tree)
- Ficus auriculata (Fig tree)
- Ficus callosa (Fig tree)
- Ficus heterophyllus (Fig tree)
- Ficus racemosa (Cluster fig tree)
- Nephelium lappaceum (Rambutan tree)
- Oroxylum indicum (Tree)
- Rubus cochinchinensis (Tree-like shrub)
- Saurauia tristyla (Tree)
- Zanthoxylum rhetsa (Prickly ash tree)
Flowering Plants (Ornamentals, Medicinal, & Wild Edibles):
- Clausena indica (Wild edible, citrus family)
- Clerodendrum chinense (Glory bower, flowering plant)
- Clerodendrum cyrtophyllum (Wild medicinal plant)
- Clerodendrum kaempferi (Ornamental plant)
- Clerodendrum philippinum (Glory flower, ornamental)
- Clerodendrum quadriloculare (Ornamental, medicinal)
- Wedelia trilobata (Creeping daisy, ground cover)
- Polygonum perfoliatum (Mile-a-minute vine, invasive)

Create species files

Create species_files a list of species and their corresponding .R files from all field visits. species_files is used for updating in the section Update based on field observations @update_model

# List of species and their corresponding `.R` files
species_files <- list(
  # === One very common species 
  bidens_pilosa = "data/effects_sd_map_bidens_pilosa.R",               
  # Bidens pilosa
  
  # === Learning in early 2024
  
  amomum_villosum = "data/effects_sd_map_amomum_villosum.R",           
  # Amomum villosum
  begonia = "data/effects_sd_map_begonia.R",                 
  # Begonia
  bischofia_javanica = "data/effects_sd_map_bischofia_javanica.R",    
  # Bischofia javanica Blum
  cheilocostus_speciosus = "data/effects_sd_map_cheilocostus_speciosus.R", 
  # Cheilocostus speciosus
  cyclosorus_acuminatus = "data/effects_sd_map_cyclosorus_acuminatus.R",
  # Cyclosorus acuminatus (Houtt)
  erythropalum_scandens = "data/effects_sd_map_erythropalum_scandens.R", 
  # Erythropalum scandens Blume
  ficus_auriculata = "data/effects_sd_map_ficus_auriculata.R", 
  # Ficus auriculata
  ficus_callosa = "data/effects_sd_map_ficus_callosa.R",     
  # Ficus callosa Willd.
  gynura_procumbens = "data/effects_sd_map_gynura_procumbens.R",        
  # Gynura procumbens (Lour.)
  melientha_suavis = "data/effects_sd_map_melientha_suavis.R",      
  # Melientha suavis
  musa_brachycarpa = "data/effects_sd_map_musa_brachycarpa.R",           
  # Musa brachycarpa Back
  ophiopogon_reptans = "data/effects_sd_map_ophiopogon_reptans.R",   
  # Ophiopogon reptans
  rubus_cochinchinensis = "data/effects_sd_map_rubus_cochinchinensis.R",     
  # Rubus cochinchinensis
  saurauia_tristyla = "data/effects_sd_map_saurauia_tristyla.R",      
  # Saurauia tristyla DC.
  schefflera_heptaphylla = "data/effects_sd_map_schefflera_heptaphylla.R", 
  # Schefflera heptaphylla
  sterculia_foetida = "data/effects_sd_map_sterculia_foetida.R",     
  # Sterculia foetida L.
  tetrapanax_papyriferus = "data/effects_sd_map_tetrapanax_papyriferus.R", 
  # Tetrapanax papyriferus
  tradescantia_pallida = "data/effects_sd_map_tradescantia_pallida.R", 
  # Tradescantia pallida
    clerodendrum_cyrtophyllum = "data/effects_sd_map_clerodendrum_cyrtophyllum.R", 
  
  # === From later visits Oct 2024
  
  # Clerodendrum cyrtophyllum
  crassocephalum_crepidioides = "data/effects_sd_map_crassocephalum_crepidioides.R", 
  # Crassocephalum crepidioides
  centella_asiatica = "data/effects_sd_map_centella_asiatica.R", 
  # Centella asiatica
  dioscorea_hamiltonii = "data/effects_sd_map_dioscorea_hamiltonii.R", 
  # Dioscorea hamiltonii
  artemisia_vulgaris = "data/effects_sd_map_artemisia_vulgaris.R", 
  # Artemisia vulgaris
  
  # === From later visits after Jan 2025
  
  gynura_sarmentosa = "data/effects_sd_map_gynura_sarmentosa.R", 
  # Gynura sarmentosa
  euodia_bodinieri = "data/effects_sd_map_euodia_bodinieri.R", 
  # Euodia bodinieri
  nephelium_lappaceum = "data/effects_sd_map_nephelium_lappaceum.R", 
  # Nephelium lappaceum
  hodgsonia_macrocarpa = "data/effects_sd_map_hodgsonia_macrocarpa.R", 
  # Hodgsonia macrocarpa
  baccaurea_ramiflora = "data/effects_sd_map_baccaurea_ramiflora.R", 
  # Baccaurea ramiflora
  cissus_hastata = "data/effects_sd_map_cissus_hastata.R", 
  # Cissus hastata
  maesa_perlarius = "data/effects_sd_map_maesa_perlarius.R", 
  # Maesa perlarius
  nephrolepis_cordifolia = "data/effects_sd_map_nephrolepis_cordifolia.R", 
  # Nephrolepis cordifolia
  piper_lolot = "data/effects_sd_map_piper_lolot.R", 
  # Piper lolot
  passiflora_foetida = "data/effects_sd_map_passiflora_foetida.R", 
  # Passiflora foetida
  plantago_major = "data/effects_sd_map_plantago_major.R", 
  # Plantago major
  zanthoxylum_rhetsa = "data/effects_sd_map_zanthoxylum_rhetsa.R", 
  # Zanthoxylum rhetsa
  costus_speciosus = "data/effects_sd_map_costus_speciosus.R", 
  # Costus speciosus
  clerodendrum_chinense = "data/effects_sd_map_clerodendrum_chinense.R", 
  # Clerodendrum chinense
  paederia_lanuginosa = "data/effects_sd_map_paederia_lanuginosa.R", 
  # Paederia lanuginosa
  clausena_indica = "data/effects_sd_map_clausena_indica.R", 
  # Clausena indica
  auricularia_heimuer = "data/effects_sd_map_auricularia_heimuer.R", 
  # Auricularia heimuer
  volvariella_volvacea = "data/effects_sd_map_volvariella_volvacea.R", 
  # Volvariella volvacea
  oroxylum_indicum = "data/effects_sd_map_oroxylum_indicum.R", 
  # Oroxylum indicum
  clerodendrum_quadriloculare = "data/effects_sd_map_clerodendrum_quadriloculare.R", 
  # Clerodendrum quadriloculare
  peperomia_pellucida = "data/effects_sd_map_peperomia_pellucida.R", 
  # Peperomia pellucida
  diplazium_esculentum = "data/effects_sd_map_diplazium_esculentum.R", 
  # Diplazium esculentum
  alpinia_xanthioides = "data/effects_sd_map_alpinia_xanthioides.R", 
  # Alpinia xanthioides
  curculigo_orchioides = "data/effects_sd_map_curculigo_orchioides.R", 
  # Curculigo orchioides
  ficus_racemosa = "data/effects_sd_map_ficus_racemosa.R", 
  # Ficus racemosa
  physalis_angulata = "data/effects_sd_map_physalis_angulata.R", 
  # Physalis angulata
  polygonum_perfoliatum = "data/effects_sd_map_polygonum_perfoliatum.R", 
  # Polygonum perfoliatum
  begonia_tonkinensis = "data/effects_sd_map_begonia_tonkinensis.R", 
  # Begonia tonkinensis
  begonia_aptera = "data/effects_sd_map_begonia_aptera.R", 
  # Begonia aptera
  canarium_bengalense = "data/effects_sd_map_canarium_bengalense.R", 
  # Canarium bengalense
  canarium_nigrum = "data/effects_sd_map_canarium_nigrum.R", 
  # Canarium nigrum
  canarium_tramvangum = "data/effects_sd_map_canarium_tramvangum.R", 
  # Canarium tramvangum
  ficus_heterophyllus = "data/effects_sd_map_ficus_heterophyllus.R", 
  # Ficus heterophyllus
  wedelia_trilobata = "data/effects_sd_map_wedelia_trilobata.R", 
  # Wedelia trilobata
  coix_lacryma_jobi = "data/effects_sd_map_coix_lacryma_jobi.R", 
  # Coix lacryma-jobi
  arundinaria_baviensis = "data/effects_sd_map_arundinaria_baviensis.R", 
  # Arundinaria baviensis
  ampelocalamus_patellaris = "data/effects_sd_map_ampelocalamus_patellaris.R", 
  # Ampelocalamus patellaris
  bambusa_blumeana = "data/effects_sd_map_bambusa_blumeana.R", 
  # Bambusa blumeana
  dendrocalamus_velutinus = "data/effects_sd_map_dendrocalamus_velutinus.R", 
  # Dendrocalamus velutinus
  momordica_cochinchinensis = "data/effects_sd_map_momordica_cochinchinensis.R", 
  # Momordica cochinchinensis
  clerodendrum_kaempferi = "data/effects_sd_map_clerodendrum_kaempferi.R", 
  # Clerodendrum kaempferi
  clerodendrum_philippinum = "data/effects_sd_map_clerodendrum_philippinum.R" 
  # Clerodendrum philippinum
)

Create `species_data`

As a last step we apply the add_species_data.R function to add the initial species data to our existing datasets (where missing, which is not usually the case).

source("functions/add_species_data.R")

Start with an empty file

# Initialize an empty dataset
species_data <- data.frame()

Apply add_species_data to add default values where any are missing.

# Loop through species and add their data
for (species in names(species_files)) {
  species_data <- add_species_data(
    existing_data = species_data,
    species_name = species,
    decision_options = decision_options,
    outcomes = outcomes,
    species_file_path = species_files[[species]]
  )
}

Load the built data and change format (for reading)

species_params <- as.data.frame(species_data)
#change to numeric
species_params$Effect <- as.numeric(species_params$Effect)
species_params$Variability_SD <- as.numeric(species_params$Variability_SD)

Clean the species level data (update botanical names).

species_params$Species <- tolower(trimws(species_params$Species))

Reverse the expenses, since higher expenses are worse, they should appear as a negative result in the plot.

Plot an aggregated plot with `aggregate_species_data`

source("functions/aggregate_species_data.R")
# Aggregate data across all species
aggregated_data <- aggregate_species_data(species_files)

# Summarize aggregated data
summary_data <- aggregated_data %>%
  group_by(Decision, Outcome) %>%
  summarize(
    MeanEffect = mean(Effect, na.rm = TRUE),
    MeanSD = mean(SD, na.rm = TRUE),
    .groups = "drop"
  )

Update based on field observations

library(ggplot2)
library(reshape2)
library(tidyverse)

Working with data from Bavi

Collected over several visits with elders and local forest food experts, chefs, teachers, youth and other villagers.

# Observed data
data <- read.csv("data/observation_data_bavi.csv")

The data is in Vietnamese so we translate it to English.

# translate names(data)

data_names <- c(
  "Family.plant", "Scientific.name", "Scientific.name.Yen.Hoai.edit",
  "Local.name", "Common.name", "Local.name...Dao.language.", "Habitat",
  "Lifeform", "Uses.To.cure.what.", "Value", "Parts used", "Usage", "Frequency
  of use" , "Number of users" , "Current reserves" , "Harvest sources" ,
  "Cultivation techniques" , "Market supply potential" , "Discovery of
  processed dishes" , "Nutritional value"
)

# Convert to snake_case (replace spaces and dots with underscores, make lowercase)
cleaned_names <- gsub("[ .]+", "_", tolower(data_names))
# Assign cleaned column names directly to the data frame
names(data) <- cleaned_names

Extract relevant information from the data

Look for all the words to be translated with unique(data$current_reserves). Add them to the dictionary in dictionaries.R and translate with translate_column.

# Source dictionary
source("functions/dictionaries.R")
# function to translate the column
source("functions/translate_column.R")

If current_reserves is Few or Not many (species is rare) we can expect positive SpeciesConservation effects with more restrictions. control_access and provide_seedlings will have a higher effect with lower SD, all the other options will likely be reduced. The effect would be similar for SustainableFoodEnv. DietaryDiversity would not be likely to change much. ParkExpenditure might also go up with more control. TraditionalKnowledge would likely not change much since it is already a rare plant. If current_reserves are Many then the effects will be the opposite.

data$current_reserves <- translate_column(data$current_reserves, translation_dict2)

Updating with parts_used data

If parts_used are roots, stems or whole plants then we can translate this as destructive and all other uses as non_destructive. The former type of use is more likely to have a negative SpeciesConservation effect with more access allowed and will shift the result.

data$parts_used <- translate_column(data$parts_used, translation_dict3)
  # destructive_harvest <- any(data$scientific_name_snake_case == species_name & data$parts_used == "destructive")

If nutritional_value is high then the effects on DietaryDiversity will be raised. If Low effects on DietaryDiversity will be lowered.

data$nutritional_value <- translate_column(data$nutritional_value, translation_dict2)

If number_of_users is Many then TraditionalKnowledge Effect would be lower and SD would be narrower for control_access. If it is Few then DietaryDiversity effects will be lower.

data$number_of_users <- translate_column(data$number_of_users, translation_dict2)

Additional functions to clean the scientific names (these are very fiddly in the rest of the steps).

  data$scientific_name <- as.character(data$scientific_name)
  
  # Function to clean and convert to snake_case
  clean_name <- function(name) {
    name <- trimws(name)  # Remove leading/trailing spaces
    name <- tolower(name)  # Convert to lowercase
    name <- gsub("[ .-]+", "_", name)  # Replace spaces and dots with underscores
    return(name)
  }
  
  # Apply function to all names
  data$scientific_name_snake_case <- sapply(data$scientific_name, clean_name, USE.NAMES = FALSE)
  
data$scientific_name_snake_case <- stringr::str_trim(data$scientific_name_snake_case)

Bayesian Updating

Use Normal-Gaussian conjugate updating for means and standard deviations. Here we use tidyverse dplyr and purrr for clean data transformations. The Effects and SDs are adjusted based on species-specific conditions. The update_effects_bayesian function assumes the prior effect is normally distributed. It incorporates new information (condition-based adjustments) using a weighted update. map() from purrr updates all the expected effects from the observations in four ways: - High/Low Nutrition: Modifies effect and SD accordingly. - User Count Impact: Adjusts effect size and uncertainty. - This Bayesian approach dynamically integrates observed data into effect estimation while preserving uncertainty, ensuring a structured update.

  # Check conditions based on the observed data

species_name = "bidens_pilosa"
rare_species <- any(data$scientific_name_snake_case == species_name &data$current_reserves == "Few")
rare_species <- any(data$scientific_name_snake_case == species_name &data$current_reserves == "Not many")
# If `current_reserves` is `Few` or `Not many` (species is rare) we can expect positive `SpeciesConservation` effects with more restrictions. `control_access` and `provide_seedlings` will have a higher effect with lower SD, all the other options will likely be reduced. The effect would be similar for `SustainableFoodEnv`. `DietaryDiversity` would not be likely to change much. `ParkExpenditure` might also go up with more control. `TraditionalKnowledge` would likely not change much  since it is already a rare plant.
common_species <- any(data$scientific_name_snake_case == species_name &data$current_reserves == "Many")
# If `current_reserves` are `Many` then the effects will be the opposite.
  destructive_harvest <- any(data$scientific_name_snake_case == species_name & data$parts_used == "destructive")
  # If `parts_used` are roots, stems or whole plants then we can translate this as `destructive` and all other uses as `non_destructive`. The former type of use is more likely to have a negative `SpeciesConservation` effect with more access allowed and will shift the result.
  high_nutrition <- any(data$scientific_name_snake_case == species_name & data$nutritional_value == "High")
  # If `nutritional_value` is high then the effects on `DietaryDiversity` will be raised.
  low_nutrition <- any(data$scientific_name_snake_case == species_name & data$nutritional_value == "Low")
  # If nutritional_value is `Low` effects on `DietaryDiversity` will be lowered.
  high_number_of_users <- any(data$scientific_name_snake_case == species_name & data$number_of_users == "Many")
  # If `number_of_users` is `Many` then `TraditionalKnowledge` Effect would be lower and SD would be narrower for `control_access`.
  few_users <- any(data$scientific_name_snake_case == species_name & data$number_of_users == "Few")
  # If number_of_users is `Few` then `DietaryDiversity` effects will be lower.

Check for species consistency

Check that the species names in the prior data species_files and the new cleaned species names in the data column scientific_name_snake_case have the same names. Check which are different and justify the missing from each list.

Species not included in first rounds

# Find species missing in prior species_files but present in observation data
setdiff(species_data_list, species_files_list)

## [1] "lactuca_indica"      "musa_acuminata"      "zingiber_officinale"
## [4] "artemisia_vulgaris_"

Missing species in species_files (present in data but not in species_files) are plants that are common but not really relevant for the question of access to forest species and so are not included in the assessment.

Species not included in later rounds

# Find extra species in prior species_files but not in observation data
setdiff(species_files_list, species_data_list)

##  [1] "amomum_villosum"        "begonia"                "cheilocostus_speciosus"
##  [4] "cyclosorus_acuminatus"  "musa_brachycarpa"       "ophiopogon_reptans"    
##  [7] "sterculia_foetida"      "tetrapanax_papyriferus" "tradescantia_pallida"  
## [10] "artemisia_vulgaris"

Extra species in species_files (present in species_files but not in data) are plants that we recognized a food use for early on in field explorations but did not learn any more about them in later visits.

Generalize the update to all species in the collection (from all trips)

Define directories

library(purrr)
library(fs)
library(stringr)

# Define directories
species_dir <- "data"
updated_dir <- "data/updated"

# Ensure updated directory exists
dir_create("data/updated")

Also assign a base directory for the updated data

# Ensure base directory
base_dir <- "data/updated"

List all the files that are in the base_dir in the new updated_species_files

# Get all updated files
updated_species_files <- list.files(base_dir, pattern = "\\.R$", full.names = TRUE)

Get all `unique_species`

Trim away trailing and leading spaces with trimws

# Get all unique species
unique_species <- unique(trimws(data$scientific_name_snake_case))  # Trim spaces

Create `process_species`

The update to the data is based on observations in the field and these were responses that need numeric interpretation. Here is how the interpretation went.

State	Effect Change	SD Change	Interpretation
Rare Species	+0.05 to +0.1	Decrease (more certainty)	Conservation actions are more effective with rare species, leading to a greater positive effect. Less uncertainty because of urgency and targeted efforts.
Common Species	-0.05 to -0.1	Increase (more uncertainty)	Conservation efforts have less impact on common species, leading to a reduced effect. Higher uncertainty due to broader variability in population responses.
Destructive Harvest	-0.05	+0.005	Overharvesting reduces conservation effectiveness. Increased uncertainty due to varied management responses.
High Nutrition	+0.05	Decrease	Higher nutrition leads to stronger dietary diversity improvements. More certainty due to known nutritional benefits.
Low Nutrition	-0.05	Increase	Poor nutrition contributes less to dietary diversity. Greater uncertainty due to variability in dietary adoption.
Many Users	-0.05	Decrease	High user numbers reduce the individual impact of knowledge sharing. More certainty in effect estimates due to larger adoption.
Few Users	-0.05	Increase	Low adoption reduces knowledge-sharing impact. Increased uncertainty due to variability in individual responses.

Create updated data for the assessment of management options for food environments and conservation.

# Function to process each species

process_species <- function(species_name) {
  
  # Construct file paths
  species_filename <- str_c("effects_sd_map_", species_name, ".R")
  species_filepath <- file.path(species_dir, species_filename)
  updated_filename <- species_filename  # Keep the same filename
  updated_filepath <- file.path(updated_dir, updated_filename)
  
  # Check if file exists
  if (!file.exists(species_filepath)) {
    message("File not found for: ", species_name)
    return(NULL)
  }
  
  message("Sourcing R script: ", species_filepath)
  
  # Source the R script to load the effect data
  source(species_filepath, local = TRUE)
  
  # Identify the loaded variable dynamically
  loaded_objects <- ls(pattern = "effects_sd_map_")
  
  if (length(loaded_objects) == 0) {
    message("No object found in: ", species_filename)
    return(NULL)
  }
  
  # Extract the loaded object dynamically
  effects_data <- get(loaded_objects[1])
  
  # Define prior knowledge-based parameters to update
  rare_species <- any(data$scientific_name_snake_case == species_name & data$current_reserves %in% c("Few", "Not many"))
  common_species <- any(data$scientific_name_snake_case == species_name & data$current_reserves == "Many")
  destructive_harvest <- any(data$scientific_name_snake_case == species_name & data$parts_used == "destructive")
  high_nutrition <- any(data$scientific_name_snake_case == species_name & data$nutritional_value == "High")
  low_nutrition <- any(data$scientific_name_snake_case == species_name & data$nutritional_value == "Low")
  high_number_of_users <- any(data$scientific_name_snake_case == species_name & data$number_of_users == "Many")
  few_users <- any(data$scientific_name_snake_case == species_name & data$number_of_users == "Few")
  
if (rare_species) {
  effects_data$SpeciesConservation$control_access$Effect <- effects_data$SpeciesConservation$control_access$Effect + 0.1
  effects_data$SpeciesConservation$control_access$SD <- pmax(0, effects_data$SpeciesConservation$control_access$SD - 0.01)

  effects_data$SpeciesConservation$provide_seedlings$Effect <- effects_data$SpeciesConservation$provide_seedlings$Effect + 0.15
  effects_data$SpeciesConservation$provide_seedlings$SD <- pmax(0, effects_data$SpeciesConservation$provide_seedlings$SD - 0.1)
}

# Adjust SpeciesConservation & SustainableFoodEnv for common species
if (common_species) {
  effects_data$SpeciesConservation$control_access$Effect <- effects_data$SpeciesConservation$control_access$Effect - 0.1
  effects_data$SpeciesConservation$control_access$SD <- pmax(0, effects_data$SpeciesConservation$control_access$SD + 0.1)

  effects_data$SpeciesConservation$provide_seedlings$Effect <- effects_data$SpeciesConservation$provide_seedlings$Effect - 0.15
  effects_data$SpeciesConservation$provide_seedlings$SD <- pmax(0, effects_data$SpeciesConservation$provide_seedlings$SD + 0.1)
}

# Adjust for destructive harvesting
if (destructive_harvest) {
  effects_data$SpeciesConservation$control_access$Effect <- effects_data$SpeciesConservation$control_access$Effect - 0.5
  effects_data$SpeciesConservation$control_access$SD <- pmax(0, effects_data$SpeciesConservation$control_access$SD + 0.5)

  effects_data$SpeciesConservation$provide_access_control_harvest$Effect <- effects_data$SpeciesConservation$provide_access_control_harvest$Effect - 0.5
  effects_data$SpeciesConservation$provide_access_control_harvest$SD <- pmax(0, effects_data$SpeciesConservation$provide_access_control_harvest$SD + 0.5)

  effects_data$SpeciesConservation$allow_access_forest_parts$Effect <- effects_data$SpeciesConservation$allow_access_forest_parts$Effect - 0.5
  effects_data$SpeciesConservation$allow_access_forest_parts$SD <- pmax(0, effects_data$SpeciesConservation$allow_access_forest_parts$SD + 0.5)

  effects_data$SpeciesConservation$allow_access_limited_time$Effect <- effects_data$SpeciesConservation$allow_access_limited_time$Effect - 0.5
  effects_data$SpeciesConservation$allow_access_limited_time$SD <- pmax(0, effects_data$SpeciesConservation$allow_access_limited_time$SD + 0.5)
}

# Adjust for high/low nutrition
if (high_nutrition) {
  effects_data$DietaryDiversity$Effect <- effects_data$DietaryDiversity$Effect + 0.05
  effects_data$DietaryDiversity$SD <- pmax(0, effects_data$DietaryDiversity$SD - 0.5)
}

if (low_nutrition) {
  effects_data$DietaryDiversity$Effect <- effects_data$DietaryDiversity$Effect - 0.5
  effects_data$DietaryDiversity$SD <- pmax(0, effects_data$DietaryDiversity$SD + 0.5)
}

# Adjust for high/few users
if (high_number_of_users) {
  effects_data$TraditionalKnowledge$control_access$Effect <- effects_data$TraditionalKnowledge$control_access$Effect - 0.5
  effects_data$TraditionalKnowledge$control_access$SD <- pmax(0, effects_data$TraditionalKnowledge$control_access$SD - 0.5)
}

if (few_users) {
  effects_data$DietaryDiversity$Effect <- effects_data$DietaryDiversity$Effect + 0.5
  effects_data$DietaryDiversity$SD <- pmax(0, effects_data$DietaryDiversity$SD - 0.5)
}


  # Assign modified data back before dumping
  assign(loaded_objects[1], effects_data)

  # Save updated data as an R script
  dump(loaded_objects[1], file = updated_filepath)
  
  message("Updated and saved: ", updated_filename)
}

# Extract species names by removing prefix and suffix
species_names <- sub("data/updated/effects_sd_map_(.*)\\.R", "\\1", updated_species_files)

Convert species_names to a named list

# Convert to a named list
updated_species_files_named <- setNames(updated_species_files, species_names)

Use walk function from purrr to apply unique species names and process_species to all species in the data.

# Apply function to all species
walk(names(updated_species_files_named), process_species)

Plot an aggregated plot with `aggregate_species_data`

Plot generate aggregated data for all species in updated_species_files_named with aggregate_species_data

### Plot an aggregated plot with `aggregate_species_data`

# source("functions/aggregate_species_data.R")
# Aggregate data across all species
# aggregated_data <- aggregate_species_data(updated_species_files_named)

library(dplyr)
library(purrr)
library(tidyr)

# # Define directory
updated_dir <- "data/updated"
#
# # List all .R files in the directory
species_files <- list.files(updated_dir, pattern = "\\.R$", full.names = TRUE)

# # Aggregate all species data
aggregated_data_updated <- map_dfr(species_files, ~ {
  source(.x, local = TRUE)  # Source the file
  species_name <- gsub(".*effects_sd_map_(.*)\\.R", "\\1", .x)  # Extract species name
  object_name <- paste0("effects_sd_map_", species_name)

  if (!exists(object_name, envir = environment())) {
    return(NULL)  # Skip if the object doesn't exist
  }

  effects_sd_map <- get(object_name, envir = environment())

  # Ensure the structure is correct
  map_dfr(names(effects_sd_map), function(outcome) {
    map_dfr(names(effects_sd_map[[outcome]]), function(decision) {
      decision_data <- effects_sd_map[[outcome]][[decision]]

      # Check if data has both Effect and SD; otherwise, skip
      if (!is.list(decision_data) || !all(c("Effect", "SD") %in% names(decision_data))) {
        return(NULL)
      }

      tibble(
        Species = species_name,
        Outcome = outcome,
        Decision = decision,
        Effect = decision_data$Effect,
        SD = decision_data$SD
      )
    })
  })
})

# Summarize aggregated data
summary_data_updated <- aggregated_data_updated %>%
  group_by(Decision, Outcome) %>%
  summarize(
    MeanEffect = mean(Effect, na.rm = TRUE),
    MeanSD = mean(SD, na.rm = TRUE),
    .groups = "drop"
  )

Apply plot_aggregate_effects to show the updated expectations across all the species in the study.

library(ggplot2)
source("functions/plot_aggregate_effects.R")
# Plot the data
plot_aggregate_effects(summary_data_updated)

The aggregated results across all species reveal distinct patterns in the effects of different decision options on key outcomes, with varying levels of confidence (represented by the standard deviation, SD). Among the interventions, provide_seedlings demonstrates the strongest overall positive effects, particularly on TraditionalKnowledge (0.43, SD = 0.05) and SustainableFoodEnv (0.28, SD = 0.04), suggesting that increasing seedling availability contributes significantly to knowledge retention and environmental sustainability. However, it also leads to a notable reduction in ParkExpenditure (-0.12, SD = 0.06), indicating potential trade-offs in resource allocation.

Other interventions, such as control_access, show mixed effects. While it improves SpeciesConservation (0.24, SD = 0.14), it appears to slightly reduce TraditionalKnowledge (-0.09, SD = 0.03), potentially due to restrictions on community engagement with the species.

Similarly, allow_access_forest_parts has a positive effect on DietaryDiversity (0.16, SD = 0.06) but has a neutral or slightly negative impact on ParkExpenditure (-0.01, SD = 0.06).

Overall, these results highlight key trade-offs between conservation, knowledge retention, and resource allocation, emphasizing the importance of balancing access policies with sustainability goals.

Forest access for sustainable food environments in Bavi

Generalized Function

Function Purpose

Simulation Workflow

Decision Options:

Estimates

Add species with simplified initial input

Decision options for sustainable forest food

Outcomes of sustainable forest food decisions

Species Classification

Create species files

Create species_data

Load the built data and change format (for reading)

Clean the species level data (update botanical names).

Plot an aggregated plot with aggregate_species_data

Update based on field observations

Working with data from Bavi

Extract relevant information from the data

Updating with parts_used data

Bayesian Updating

Check for species consistency

Species not included in first rounds

Species not included in later rounds

Generalize the update to all species in the collection (from all trips)

Define directories

Get all unique_species

Create process_species

Plot an aggregated plot with aggregate_species_data

Create `species_data`

Plot an aggregated plot with `aggregate_species_data`

Get all `unique_species`

Create `process_species`

Plot an aggregated plot with `aggregate_species_data`